Multiple Linear Regression¶
Data Preprocessing¶
It is a data mining technique that transforms raw data into an understandable format. Raw data(real world data) is always incomplete and that data cannot be sent through a model. That would cause certain errors. That is why we need to preprocess data before sending through a model.¶
Steps in Data Preprocessing These are the steps:
- Import libraries
- Import dataset
- Finding for missing values
- Encoding categorical data
- Data splitting
- Feature Scaling
1. Importing Libraries¶
In [1]:
#Importing the Libraries
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
import statsmodels.api as sm
from sklearn.linear_model import LinearRegression
from sklearn.datasets import load_boston
from sklearn.model_selection import train_test_split
from statsmodels.stats.outliers_influence import variance_inflation_factor
2. Importing Dataset¶
In [2]:
# Gather Data
boston_dataset = load_boston()
In [3]:
type(boston_dataset)
Out[3]:
In [4]:
# boston_dataset
dir(boston_dataset)
Out[4]:
In [5]:
print(boston_dataset.DESCR)
In [6]:
print(boston_dataset.data)
In [7]:
print(boston_dataset.feature_names)
In [8]:
print(boston_dataset.filename)
In [9]:
print(boston_dataset.target)
In [10]:
data = pd.DataFrame(data = boston_dataset.data, columns = boston_dataset.feature_names)
data.shape
Out[10]:
In [11]:
data.head()
Out[11]:
In [12]:
data['PRICE'] = boston_dataset.target
In [13]:
data.head()
Out[13]:
In [14]:
data.tail()
Out[14]:
In [15]:
data.count()
Out[15]:
In [16]:
data.info()
3. Checking missing value¶
In [17]:
data.isnull()
Out[17]:
In [18]:
data.isnull().any()
Out[18]:
In [19]:
data.isnull().sum()
Out[19]:
5. Data splitting¶
In [20]:
prices = data['PRICE']
features = data.drop('PRICE', axis = 1)
In [21]:
prices
Out[21]:
In [22]:
features
Out[22]:
In [23]:
# train, test split......
X_train, X_test, y_train, y_test = train_test_split(features, prices, test_size=0.2, random_state = 10)
In [24]:
X_train
Out[24]:
In [25]:
y_train
Out[25]:
In [26]:
X_test
Out[26]:
In [27]:
y_test
Out[27]:
In [28]:
print("X_train Size :",len(X_train))
print("Y_train Size :",len(y_train))
print("X_test Size :",len(X_test))
print("Y_test Size :",len(y_test))
print("Train Size :", (len(X_train)/len(features))*100)
print("Train Size :", (len(X_test)/len(features))*100)
Data Visuvalization¶
In [29]:
plt.hist(data['PRICE'])
plt.show()
In [30]:
plt.hist(data['PRICE'], bins = 30)
plt.show()
In [31]:
plt.hist(data['PRICE'], bins = 30, ec ='black', color = '#800000', alpha = 0.5)
plt.xlabel('Price in 1000\'s')
plt.ylabel('Nr of Houses')
plt.show()
In [32]:
sns.distplot(data['PRICE'])
plt.show()
In [33]:
plt.figure(figsize= (10,6))
sns.distplot(data['PRICE'], color = 'red', bins = 50)
plt.show()
In [34]:
plt.figure(figsize= (10,6))
sns.distplot(data['PRICE'], color = 'red', bins = 50, hist = False)
plt.show()
In [35]:
plt.figure(figsize= (10,6))
sns.distplot(data['PRICE'], color = 'red', bins = 50, kde = False)
plt.show()
In [36]:
plt.figure(figsize= (10,6))
plt.hist(data['RM'], ec = 'black', color = 'pink')# removed bins
plt.xlabel('Average Number of Rooms')
plt.ylabel('Nr of Houses')
plt.show()
In [37]:
data['RM'].mean()
Out[37]:
In [38]:
plt.figure(figsize= (10,6))
plt.hist(data['RAD'], ec = 'black', color = 'green', alpha = 0.5)
plt.xlabel('Accessibility to Radial Highways')
plt.ylabel('Nr of Houses')
plt.show()
In [39]:
data['RAD'].value_counts()
Out[39]:
In [40]:
accessibility = data['RAD'].value_counts()
type(accessibility)
Out[40]:
In [41]:
print(accessibility.index)
In [42]:
print(accessibility.values)
In [43]:
# my bar chart
plt.figure(figsize= (10,6))
plt.bar(accessibility.index, accessibility )
plt.show()
In [44]:
plt.figure(figsize= (10,6))
plt.bar(accessibility.index, accessibility, color ='brown', ec = 'blue', width = 0.5)# changing the width, color and ec
plt.xlabel('Accessibility index to Highways')
plt.ylabel('Nr of Houses')
plt.show()
In [45]:
data['CHAS'].value_counts()
Out[45]:
In [46]:
plt.figure(figsize= (10,6))
plt.hist(data['CHAS'], ec = 'black', color = 'green', alpha = 0.5)
plt.xlabel('Near Charles River')
plt.ylabel('Nr of Houses')
plt.show()
In [47]:
# Smallest value, Largest Value, Mean Value, Median Value..
data['PRICE'].min()
Out[47]:
In [48]:
data['PRICE'].max()
Out[48]:
In [49]:
data['PRICE'].mean()
Out[49]:
In [50]:
data['PRICE'].median()
Out[50]:
In [51]:
data.min()
Out[51]:
In [52]:
data.max()
Out[52]:
In [53]:
data.mean()
Out[53]:
In [54]:
data.median()
Out[54]:
In [55]:
data.describe()
Out[55]:
In [56]:
data['RM'].describe()
Out[56]:
In [57]:
data['PRICE'].corr(data['RM'])
Out[57]:
In [58]:
# Price and PTRATIO pupil-teacher ratio by town
# pupil-teacher ratio ???
# pupil - teacher ratio
# pupil/teacher ??? .....21/1........pupil ratio 21
# smaller the pupil - teacher ratio the better?
data['PRICE'].corr(data['PTRATIO'])
Out[58]:
In [59]:
data.corr()
Out[59]:
In [60]:
type(data.corr())
Out[60]:
In [61]:
masker = np.zeros_like(data.corr())
In [62]:
masker
Out[62]:
In [63]:
triangle_indices = np.triu_indices_from(masker, k = -2)
masker[triangle_indices] = True
masker
Out[63]:
In [64]:
# Heat map
plt.figure(figsize =(16,10))
sns.heatmap(data.corr())
plt.show()
In [65]:
plt.figure(figsize =(10,8))# reduced the size
sns.heatmap(data.corr())
plt.xticks(fontsize = 10)#adjusted the fonts on the x axis
plt.yticks(fontsize = 10)# adjusted the font on the y axis
plt.show()
In [66]:
plt.figure(figsize =(10,8))
sns.heatmap(data.corr(), mask = masker)
plt.xticks(fontsize = 10)
plt.yticks(fontsize = 10)
plt.show()
In [67]:
plt.figure(figsize =(10,8))
sns.heatmap(data.corr(), mask = masker, annot=True)
plt.xticks(fontsize = 10)
plt.yticks(fontsize = 10)
plt.show()
In [68]:
plt.figure(figsize =(16,10))
sns.heatmap(data.corr(), mask = masker, annot=True,annot_kws={'size': 14} )
plt.xticks(fontsize = 10)
plt.yticks(fontsize = 10)
plt.show()
In [69]:
# Advanced Visualisations
nox_dis_corr = round(data['NOX'].corr(data['DIS']), 3)
nox_dis_corr
plt.figure(figsize = (16, 8))
plt.scatter(x=data['DIS'], y =data['NOX'], alpha = 0.5, color = 'green', s= 80)
plt.title(f' Distance from Empy centres versus Pollution, Corr({nox_dis_corr})', fontsize = 25, color ='brown' )
plt.xlabel( 'Distance from Employment Centres...DIS', fontsize =14, color = 'black' )
plt.ylabel( 'Level of Pollutants Nitric Oxide...NOX', fontsize =14, color = 'black' )
plt.show()
In [70]:
price_rad_corr = round(data['PRICE'].corr(data['RAD']), 3)
price_rad_corr
plt.figure(figsize = (16, 8))
plt.scatter(x=data['PRICE'], y =data['RAD'], alpha = 0.5, color = 'green', s= 80)
# plt.title(f' Distance from Empy centres versus Pollution, Corr({nox_dis_corr})', fontsize = 25, color ='brown' )
# plt.xlabel( 'Distance from Employment Centres...DIS', fontsize =14, color = 'black' )
# plt.ylabel( 'Level of Pollutants Nitric Oxide...NOX', fontsize =14, color = 'black' )
plt.show()
In [71]:
price_rad_corr
Out[71]:
In [72]:
rad_price_corr = round(data['RAD'].corr(data['PRICE']), 3)
rad_price_corr
plt.figure(figsize = (16, 8))
plt.scatter(x=data['RAD'], y =data['PRICE'], alpha = 0.5, color = 'green', s= 80)
# plt.title(f' Distance from Empy centres versus Pollution, Corr({nox_dis_corr})', fontsize = 25, color ='brown' )
# plt.xlabel( 'Distance from Employment Centres...DIS', fontsize =14, color = 'black' )
# plt.ylabel( 'Level of Pollutants Nitric Oxide...NOX', fontsize =14, color = 'black' )
plt.show()
In [73]:
# using seaborn to plot
sns.set() # resets any previous styling
sns.set_style('dark')
sns.jointplot(x = data['DIS'], y = data['NOX'], color = 'red', joint_kws = {'alpha':0.2})
plt.show()
In [74]:
sns.set() # to reset any previous styling
sns.set_style('whitegrid')# to set style.
plt.figure(figsize = (16, 8))
sns.scatterplot(x=data['DIS'], y =data['NOX'], color = 'red', alpha = 0.5) # joint plot creates scatter plot
plt.show()
In [75]:
sns.set()
sns.set_style('dark')
sns.jointplot(x=data['DIS'], y =data['NOX'], kind = 'hex' ) # kind is hex, darker hexagons in regions of higher density
plt.show()
In [76]:
# TAX and RAD
sns.set() # to reset any previous styling
sns.set_style('dark')# to set style.
sns.jointplot(x=data['TAX'], y =data['RAD'], color = 'darkred', joint_kws = {'alpha' : 0.5} ) # jointplot gives scatter plot
plt.show()
In [77]:
# regression plot...lmplot....
sns.lmplot(x ='TAX', y = 'RAD', data = data)
plt.show()
In [78]:
sns.lmplot(x ='TAX', y = 'RAD', data = data, height = 7)
plt.show()
In [79]:
sns.lmplot(x ='RM', y = 'PRICE', data = data, height =10, aspect = 2, ci= None )
plt.show()
In [80]:
sns.pairplot(data)
plt.show()
In [81]:
sns.pairplot(data, kind = 'reg' , plot_kws = {'line_kws':{'color': 'red'}})
plt.show()
In [82]:
regr = LinearRegression()
regr.fit(X_train, y_train) # I am using the training data to train the model
Out[82]:
In [83]:
regr.coef_
Out[83]:
In [84]:
regr.intercept_
Out[84]:
In [85]:
pd.DataFrame(data = regr.coef_, index = X_train.columns, columns = ['Coef'])
Out[85]:
In [86]:
regr.score(X_train, y_train) # score on train data, that the model has seen
Out[86]:
In [87]:
regr.score(X_test, y_test) # score on the test data the model has not seen before
Out[87]:
In [88]:
plt.figure(figsize= (16,8))
plt.hist(data['PRICE'], bins =50, ec = 'black', color = 'blue', alpha = 0.7)# bins 50...setting color by hexcode
plt.xlabel('Price in 1000\'s')
plt.ylabel('Nr of Houses')
plt.show()
In [89]:
# The plot above has a skew.....Its not balanced by the middle value
# Lets find the skew..
data['PRICE'].skew()
Out[89]:
In [90]:
#skew = 0 but skew = 1.1
In [91]:
y_log = np.log(data['PRICE'])
y_log.head()
Out[91]:
In [92]:
y_log.tail()
Out[92]:
In [93]:
data['PRICE'].skew()
Out[93]:
In [94]:
y_log.skew()
Out[94]:
In [95]:
sns.set()
plt.figure(figsize=(20,6))
sns.distplot(y_log)
plt.show()
In [96]:
sns.set()
plt.figure(figsize=(20,6))
sns.distplot(y_log)
plt.title(f'Log price: The skew is {round(y_log.skew(),2)}', fontsize = 25)# note the f string used to print the y_log.skew()
plt.xlabel('LOG PRICES ', fontsize =14, color = 'green' )
plt.show()
In [97]:
#plot before log transformation for comparison
skew=round(data['PRICE'].skew(),2)
sns.set()
plt.figure(figsize=(20,6))
sns.distplot(data['PRICE'])
plt.title(f'Actual price: The skew is {skew}', fontsize = 25)
plt.xlabel( 'ACTUAL PRICES ', fontsize =14, color = 'green' )
plt.show()
In [98]:
# doing a scatter plot with LSTAT vs PRICE .....on untransformed data...using sns......
sns.lmplot(x='LSTAT', y='PRICE', data=data, height=10, aspect =2,
scatter_kws={'alpha': 0.6}, line_kws={'color':'darkred'})
plt.show()
In [99]:
# doing a scatter plot with LSTAT vs PRICE .....on log transformed data...using sns......
# transforming the price into log price
data
Out[99]:
In [100]:
transformed_data = features
transformed_data['LOG_PRICE'] = y_log
transformed_data
Out[100]:
In [101]:
sns.lmplot(x='LSTAT', y='LOG_PRICE', data=transformed_data, height=10, aspect =2,
scatter_kws={'alpha': 0.6}, line_kws={'color':'red'})
plt.show()
In [102]:
prices = np.log(data['PRICE'])
features = data.drop('PRICE', axis=1)
X_train, X_test, y_train, y_test = train_test_split(features, prices,
test_size=0.2, random_state=10)
regr = LinearRegression()
regr.fit(X_train, y_train)
Out[102]:
In [103]:
LinearRegression(copy_X=True, fit_intercept=True, n_jobs=None, normalize=False)
Out[103]:
In [104]:
print('Training data r-squared:', regr.score(X_train, y_train))
In [105]:
print('Test data r-squared:', regr.score(X_test, y_test))
In [106]:
pd.DataFrame(data=regr.coef_, index=X_train.columns, columns=['coef'])
Out[106]:
In [107]:
print('Intercept', regr.intercept_)# intercept after log transformation of prices
In [108]:
# CHAS 0.080331.....(got after the prices were log transformed)
# how to retranform?
np.e**0.080331
Out[108]:
In [109]:
# using the stats model to get the p value...
results = sm.OLS(y_train, sm.add_constant(X_train)).fit()
results.params
Out[109]:
In [110]:
results.pvalues
Out[110]:
In [111]:
round(results.pvalues,3)
Out[111]:
In [112]:
# from the above INDUS and AGE have p values greater than 0.05
# Multicollinearity....VIF......Variance inflation factor.....
In [113]:
variance_inflation_factor(exog=sm.add_constant(X_train).values, exog_idx=10)
Out[113]:
In [114]:
variance_inflation_factor(exog=sm.add_constant(X_train).values, exog_idx=13)
Out[114]:
In [115]:
# write a loop for gettig all the VIF...and if any feature has a VIF of > 10.....then I have a pblm of multi collinearity
X_train
Out[115]:
In [116]:
X_incl_const = sm.add_constant(X_train)
X_incl_const
Out[116]:
In [117]:
# write a for loop to get the VIF of all features .....Home work
variance_inflation_factor(exog=sm.add_constant(X_train).values, exog_idx=13)
Out[117]:
In [118]:
X_incl_const.shape[1]
Out[118]:
In [119]:
vif = [] # empty list
for i in range (X_incl_const.shape[1]):
vif.append(variance_inflation_factor(exog=sm.add_constant(X_train).values, exog_idx=i))
print(vif)
In [120]:
pd.DataFrame({'coef_name': X_incl_const.columns, 'vif': np.around(vif, 2)})
Out[120]:
In [121]:
# BIC Test
# This model is using all the features and the log prices...
In [122]:
# model_1
X_incl_const = sm.add_constant(X_train)
model_1 = sm.OLS(y_train, X_incl_const)
results_1 = model_1.fit()
results_1.params
Out[122]:
In [123]:
round(results_1.pvalues,3)
Out[123]:
In [124]:
round(results_1.bic,2)
Out[124]:
In [125]:
# Another model_2 after dropping the INDUS(high p value) and using the log prices..
X_incl_const = sm.add_constant(X_train)
X_incl_const = X_incl_const.drop(['INDUS'], axis = 1)
model_2 = sm.OLS(y_train, X_incl_const)
results_2 = model_2.fit()
results_2.params
Out[125]:
In [126]:
round(results_2.pvalues,3)
Out[126]:
In [127]:
round(results_2.bic,2)
Out[127]:
In [128]:
# third model_3 after dropping the INDUS and the AGE, and using the log prices..
In [129]:
X_incl_const = sm.add_constant(X_train)
X_incl_const = X_incl_const.drop(['INDUS', 'AGE'],axis = 1)
model_3 = sm.OLS(y_train, X_incl_const)
results_3 = model_3.fit()
results_3.params
Out[129]:
In [130]:
round(results_3.pvalues,3)
Out[130]:
In [131]:
round(results_3.bic,2)
Out[131]:
In [132]:
# Home work is....from the models model_1, model_2, model_3 .....get the results and make a DataFrame like above..
X_incl_const = sm.add_constant(X_train)
model_1 = sm.OLS(y_train, X_incl_const)
results_1 = model_1.fit()
results_1.params
results_1.pvalues
model_1_df = pd.DataFrame({'M1_Coef': results_1.params, 'M1_Pvalues' : round(results_1.pvalues,3)})
model_1_df
Out[132]:
In [133]:
X_incl_const = sm.add_constant(X_train)
model_2 = sm.OLS(y_train, X_incl_const.drop(['INDUS'], axis =1))
results_2 = model_2.fit()
results_2.params
results_1.pvalues
model_2_df = pd.DataFrame({'M2_Coef': results_2.params, 'M2_Pvalues' : round(results_2.pvalues,3)})
model_2_df
Out[133]:
In [134]:
X_incl_const = sm.add_constant(X_train)
model_3 = sm.OLS(y_train, X_incl_const.drop(['INDUS', "AGE"], axis =1))
results_3 = model_3.fit()
results_3.params
results_3.pvalues
model_3_df = pd.DataFrame({'M3_Coef': results_3.params, 'M3_Pvalues' : round(results_3.pvalues,3)})
model_3_df
Out[134]:
In [135]:
frames = [model_1_df, model_2_df,model_3_df]
pd.concat(frames, axis =1)
Out[135]:
In [136]:
print(f'BIC WITH ALL:: {results_1.bic}: W/O INDUS {results_2.bic}: W/O INDUS AND AGE {results_3.bic} ')
In [137]:
print(f'RSQD:: WITH ALL {(results_1.rsquared)}: W/O INDUS {results_2.rsquared}: W/O INDUS,AGE {results_3.rsquared}')
In [138]:
prices = np.log(data['PRICE']) # Use log prices
features = data.drop(['PRICE', 'INDUS', 'AGE'], axis=1)# dropping the INDUS , AGE also along with PRICE which has been moved
X_train, X_test, y_train, y_test = train_test_split(features, prices,
test_size=0.2, random_state=10)
# Using Statsmodel
results = sm.OLS(y_train, sm.add_constant(X_train)).fit()# obtaining the regression coeffs including the constant(intercept).
# commented out lines below can be used to look at the residuals.
residuals = y_train - results.fittedvalues
results.resid
Out[138]:
In [139]:
# Graph of Actual vs. Predicted Prices
corr = round(y_train.corr(results.fittedvalues), 2)# correlation rounded
plt.figure(figsize=(20,6))
plt.scatter(x=y_train, y=results.fittedvalues, c='navy', alpha=0.6)# gives a scatter plot
plt.plot(y_train, y_train, color='red')#the "ideal regression" line......not a reality
plt.xlabel(' log prices $y _i$', fontsize=20)
plt.ylabel('Prediced log prices $\hat y _i$', fontsize=20)
plt.title(f'Log prices vs Predicted log prices: $y _i$ vs $\hat y_i$ (Corr {corr})', fontsize=20)
plt.show()
In [140]:
# Actual Prices vs Predicted Prices in 1000's of Dollars
plt.figure(figsize=(20,6))
plt.scatter(x=np.e**y_train, y=np.e**results.fittedvalues, c='blue', alpha=0.6)
plt.plot(np.e**y_train, np.e**y_train, color='brown')
plt.xlabel('Actual prices in Dollars $y _i$', fontsize=20)
plt.ylabel('Prediced prices in Dollars $\hat y _i$', fontsize=20)
plt.title(f'Actual vs Predicted prices: $y _i$ vs $\hat y_i$ (Corr {corr})', fontsize=20)
plt.show()
In [141]:
# Residuals vs Predicted values
plt.figure(figsize=(20,6))
plt.scatter(x=results.fittedvalues, y=results.resid, c='navy', alpha=0.6)
plt.xlabel('Predicted log prices $\hat y _i$', fontsize=20)
plt.ylabel('Residuals', fontsize=20)
plt.title('Residuals vs Fitted Values', fontsize=20)
plt.show()
# Mean Squared Error & R-Squared
reduced_log_mse = round(results.mse_resid, 3)# for later use and comparision
reduced_log_rsquared = round(results.rsquared, 3) # for later use and comparision
In [142]:
# Distribution of Residuals (log prices) - checking for normality
resid_mean = round(results.resid.mean(), 3)
resid_skew = round(results.resid.skew(), 3)
In [143]:
print(resid_mean, '\n');
print(resid_skew) # printing the mean and skew, both a close to zero
In [144]:
plt.figure(figsize= (20,6))
sns.distplot(results.resid, color='navy')# using the distplot of seaborn
plt.title(f'Log price model: residuals Skew ({resid_skew}) Mean ({resid_mean})')
plt.show()
In [145]:
# Using the original model with all the features and normal prices generate:
# Plot of actual vs predicted prices (incl. correlation)
# Plot of residuals vs. predicted prices
# Plot of distribution of residuals (incl. skew)
# Analysis of results.
In [146]:
# Original model: normal prices & all features
prices = data['PRICE']
features = data.drop(['PRICE'], axis=1)
X_train, X_test, y_train, y_test = train_test_split(features, prices,
test_size=0.2, random_state=10)
X_incl_const = sm.add_constant(X_train)
model = sm.OLS(y_train, X_incl_const)
results = model.fit()
# Graph of Actual vs. Predicted Prices
corr = round(y_train.corr(results.fittedvalues), 2)
plt.figure(figsize=(20,6))
plt.scatter(x=y_train, y=results.fittedvalues, c='indigo', alpha=0.6)
plt.plot(y_train, y_train, color='cyan')
plt.xlabel('Actual prices 000s $y _i$', fontsize=14)
plt.ylabel('Prediced prices 000s $\hat y _i$', fontsize=14)
plt.title(f'Actual vs Predicted prices: $y _i$ vs $\hat y_i$ (Corr {corr})', fontsize=17)
plt.show()
# Residuals vs Predicted values
plt.figure(figsize=(20,6))
plt.scatter(x=results.fittedvalues, y=results.resid, c='indigo', alpha=0.6)
plt.xlabel('Predicted prices $\hat y _i$', fontsize=14)
plt.ylabel('Residuals', fontsize=14)
plt.title('Residuals vs Fitted Values', fontsize=17)
plt.show()
# Residual Distribution Chart
plt.figure(figsize=(20,6))
resid_mean = round(results.resid.mean(), 3)
resid_skew = round(results.resid.skew(), 3)
sns.distplot(results.resid, color='indigo')
plt.title(f'Residuals Skew ({resid_skew}) Mean ({resid_mean})')
plt.show()
# Mean Squared Error & R-Squared
full_normal_mse = round(results.mse_resid, 3)# for later use
full_normal_rsquared = round(results.rsquared, 3) # for later use
In [147]:
# Model Omitting Key Features using log prices
prices = np.log(data['PRICE'])
features = data.drop(['PRICE', 'INDUS', 'AGE', 'LSTAT', 'RM', 'NOX', 'CRIM'], axis=1)
X_train, X_test, y_train, y_test = train_test_split(features, prices,
test_size=0.2, random_state=10)
X_incl_const = sm.add_constant(X_train)
model = sm.OLS(y_train, X_incl_const)
results = model.fit()
# Graph of Actual vs. Predicted Prices
corr = round(y_train.corr(results.fittedvalues), 2)
plt.scatter(x=y_train, y=results.fittedvalues, c='#e74c3c', alpha=0.6)
plt.plot(y_train, y_train, color='cyan')
plt.xlabel('Actual log prices $y _i$', fontsize=14)
plt.ylabel('Predicted log prices $\hat y _i$', fontsize=14)
plt.title(f'Actual vs Predicted prices with omitted variables: $y _i$ vs $\hat y_i$ (Corr {corr})', fontsize=17)
plt.show()
# Residuals vs Predicted values
plt.scatter(x=results.fittedvalues, y=results.resid, c='#e74c3c', alpha=0.6)
plt.xlabel('Predicted prices $\hat y _i$', fontsize=14)
plt.ylabel('Residuals', fontsize=14)
plt.title('Residuals vs Fitted Values', fontsize=17)
plt.show()
# Mean Squared Error & R-Squared
omitted_var_mse = round(results.mse_resid, 3)# for later use
omitted_var_rsquared = round(results.rsquared, 3)# for later use
In [148]:
pd.DataFrame({'R-Squared': [reduced_log_rsquared, full_normal_rsquared, omitted_var_rsquared],
'MSE': [reduced_log_mse, full_normal_mse, omitted_var_mse]},
index=['Reduced Log Model', 'Full Normal Price Model', 'Omitted Var Model'])
Out[148]:
In [149]:
pd.DataFrame({'R-Squared': [reduced_log_rsquared, full_normal_rsquared, omitted_var_rsquared],
'MSE': [reduced_log_mse, full_normal_mse, omitted_var_mse],
'RMSE': np.sqrt([reduced_log_mse, full_normal_mse, omitted_var_mse])},# taking the square root of MSE to get RMSE.
index=['Reduced Log Model', 'Full Normal Price Model', 'Omitted Var Model'])
Out[149]:
In [150]:
# use reduced log model...
# Estimated price is 30,000 DOllars for a house...
# calculate the upper bound and lower bound of the prices..
# from the above DF MSE = 0.035
In [151]:
np.sqrt(reduced_log_mse)
Out[151]:
In [152]:
2*np.sqrt(reduced_log_mse)
Out[152]:
In [153]:
upper_bound = np.log(30) + 2*np.sqrt(reduced_log_mse)
In [154]:
(np.e**upper_bound)*1000
Out[154]:
In [155]:
lower_bound = np.log(30) - 2*np.sqrt(reduced_log_mse)
In [156]:
(np.e**lower_bound)*1000
Out[156]:
In [ ]: